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ABSTRACT 

Viscous Keplerian discs become sub-Keplerian close to a black hole since they pass 
through sonic points before entering into it. We study the time evolution of polytropic 
viscous accretion discs (both in one and two dimensional flows) using Smoothed Parti¬ 
cle Hydrodynamics. We discover that for a large region of the parameter space, when 
the flow viscosity parameter is less than a critical value, standing shock waves are 
formed. If the viscosity is very high then the shock disappears. In the intermediate 
viscosity the disc oscillates very significantly in viscous time-scale. Our simulations 
indicate that these centrifugally supported high density region close to a black hole 
plays an active role in the flow dynamics, and consequently, the radiation dynamics. 

Key words: accretion, accretion discs - black hole physics - shock-waves - hydro¬ 
dynamics 


I. INTRODUCTION 

One of the convincing ways to understand the nature of 
accretion processes on black holes is to use numerical tech¬ 
niques. This is because the steady solutions of analytical 
models are too simplistic and they usually concentrate on 
steady state solutions. The works present in the literature, 
which study the possibility of sub-Keplerian accretion flows 
are mainly carried out for inviscid flows (Hawley, Smarr 
& Wilson, 1984; Molteni, Lanzafame & Chakrabarti, 1994; 
[hereafter MLC94] Nobuto & Hanawa, 1994; Molteni, Ryu 
& Chakrabarti, 1996 [hereafter MRC96]; Ryu, Chakrabarti 
& Molteni, 1997 [hereafter RCM97[) although some simu¬ 
lations of viscous Keplerian discs are also present (Eggum, 
Coroniti & Katz, 1988). In Chakrabarti & Molteni (1995, 
hereafter CM95; and references therein), the time evolu¬ 
tion of one dimensional isothermal viscous transonic flows 
(VTFs) has been studied. It was particularly noted that 
in some region of the parameter space, when the viscos¬ 
ity is smaller, the flow developed centrifugal pressure sup¬ 
ported shock waves as is expected from the steady solutions 
(Chakrabarti, 1990; hereafter C90). As viscosity is increased, 
the angular momentum transport rate in the post-shock re¬ 
gion is enhanced compared to the rate in the pre-shock flow 
and as a result the shock propagates outwards, and the flow 
becomes sub-sonic. Further out, the flow can join with a 
Keplerian accretion disc. 

In the present paper, we study the behaviour of the 
sub-Keplerian VTFs when the flow is neither isothermal nor 


restricted only to the equatorial plane as in CM95. We ex¬ 
amine the behaviour of the more general polytropic thick, 
viscous accretion discs. We And that even in two dimensional 
thick discs, shocks form and the steady shock location in¬ 
creases with viscosity as in one dimensional study of CM95. 
We also observe, to our surprise, that beyond a critical vis¬ 
cosity, when the steady shock is not expected, the flow forms 
an unsteady shock which periodically evacuates the disc. 
Beyond another critical viscosity (keeping other parameters 
unchanged), where the flow passes only through the inner 
sonic point, the shock disappears and only the smooth sub- 
Keplerian disc (originating from a Keplerian disc) remains. 
The presence oscillating solutions are clear indications of 
transitions from one topology of solutions to another. Pres¬ 
ence of both steady and oscillating solutions is significant 
in view of the fact that spectral observations of black hole 
candidates do show quasi-periodic oscillations in between 
various stationary states. Simply put, two critical viscosities 
in the above context come into being for the following rea¬ 
sons: Beyond a critical viscosity, the shock conditions cease 
to be fulfilled even when the flow has two saddle type sonic 
points. The flow passes through both, using unsteady shocks. 
Beyond yet another, and higher value of critical viscosity, 
the outer saddle type sonic point cease to exist altogether 
and flow smoothly passes through the inner sonic point. The 
general behaviour of inviscid and viscous accretion flows and 
how these critical viscosities can be identified, are discussed 
in C90 and Chakrabarti (1996, hereafter C96) respectively. 
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The organization of the present paper is the following: 
In the next Section, we briefly discuss the theory of viscous 
transonic discs. In §3, we present the simulation results in 
one dimensional viscous discs. In §4, we present the results in 
two dimensional viscous discs. In §5, we discuss astrophysical 
importance of our solutions and present concluding remarks. 

2. MODEL EQUATIONS 

We model the nonrotating central compact object us¬ 
ing Paczyhski & Wiita (1980) potential. The inclusion of 
viscosity in the form of Shakura-Sunyaev prescription (or, 
other viscosity prescriptions as in Chakrabarti & Molteni, 
1995) implies that the specific angular momentum A is not 
constant everywhere, rather it is efficiently transported out¬ 
ward at the rate determined by the magnitude of viscosity. 

We measure all distances, velocities and timescales in 
units of 2GM/c^, c and 2GM/c? respectively. Below, we 
provide the Lagrangean formulae for the two-dimensional 
fluid dynamics equations for SPH in cylindrical coordinates. 
The mass conservation equation is: 
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The radial momentum equation becomes. 
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The vertical momentum equation becomes. 
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For the Shakura Sunyaev turbulent viscosity we have used, 
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where the vertical thickness Tidisc. is estimated from the as¬ 
sumption of the vertical equilibrium condition: 
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One may use the conventional energy equation as, 
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where, e is the specihc thermal energy, and <I> is the usual 
viscous dissipation term. However, it turns out that a better 
accuracy is achieved if the sum total of specihc kinetic and 
thermal energies are used instead. This is because this quan¬ 
tity is exactly conserved (since it can be put in a symmetric 
form so that pair of particles exchange equal amounts of en¬ 
ergy, e.g. Monaghan 1985). So we prefer to use the following 
version of the energy equation: 



with 


P„ 

-Vv + v- 


P 



+ -V 


P 




-ivP+g 

P 


where P = (7 — 1) pe is the equation of state of ideal gas, 
g is the gravitational acceleration, r : v is the vector re¬ 
sulting from the contraction of the stress tensor with the 
velocity vector. We include only Trp (namely, the ref) com¬ 
ponent) since it is the dominant contributor to the vis¬ 
cous stress. We choose for simplicity the pseudo-Newtonian 
potential proposed by Paczyhski & Wiita (1980), where 
<I>(r, 9) = l/2(r — 1). A complete steady solution requires the 
equations of energy, angular momentum and mass conserva¬ 
tion supplied by transonic conditions at the critical points 
and the Rankine-Hugoniot conditions at the shock. 

Since a black hole accretion is necessarily supersonic 
at the horizon (Chakrabarti, 1990), and since the flows are 
sub-Keplerian at the sonic point even for generalized accre¬ 
tion models (C96), a flow must deviate from a Keplerian 
disc much before it enters into a hole. We verify this by per¬ 
forming numerous numerical simulations. If the viscosity is 
smaller than a critical value, and if the shock conditions are 
satisfied, then the shocks form, otherwise the shock becomes 
unstable and propagate outwards to make a subsonic disc 
which joins with a Keplerian disc. The reason for the shock 
to propagate outwards is that the rate of angular momentum 
transport in the pre and post-shock flows become completely 
different. In particular, the transport rate is higher in the 
subsonic post-shock region. Even when the shock is absent, 
the centrifugal barrier supported region close to the black 
hole behave very similar to the post-shock region, and hence 
the emission properties, of this shock-free disc also becomes 
similar (Chakrabarti, 1997). 

In this context, the difference between the approach 
to obtain the steady solutions (as in C90, or in C96) and 
time dependent solutions is worth noting. While obtaining 
a steady solution, since the flow must pass through a sonic 
point, it is customary to choose it’s location and the spe¬ 
cific angular momentum of the flow at this point. The rest 
of the properties of the flow (such as the location where the 
flow deviates from a Keplerian disc, it’s injection velocities 
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etc.) are determined as eigenvalues. In the time-dependent 
studies, on the other hand, we hx the injected parameters 
(such as the specific energy and angular momentum) while 
the location of the sonic point and the angular momentum 
at the sonic point are determined automatically as eigen¬ 
values. That these two approaches produce the same result 
(whenever steady state solutions exist) has been verified re¬ 
peatedly (See, first paragraph of the Introduction for refer¬ 
ence.). In the present paper also, we shall fix the parameters 
of the injecting matter at the outer boundary. 

The numerical simulations we perform are done using 
Smoothed Particle Hydrodynamics (SPH). The code has 
been tested several times against theoretical works (e.g., 
CM93, CM95; Sponholz & Molteni, 1994; MRC96; RCM97) 
and is found to be sufficiently accurate as far as the angular 
momentum conservation is concerned. In the next section 
we present results in one dimensions and in §4, we present 
results in 2 dimensional disc. 

2. SIMULATION RESULTS IN ONE DIMENSIONS 

Figure 1 shows stationary shock locations as a function 
of Shakura-Sunyaev viscosity parameter a in thin, rotating, 
accreting flows. The viscosity coefficient u is related to the 
viscosity parameter a by (see, Shakura & Sunyaev, 1973) 
fi = apaZdisc- In upper left panel, we plot density p{r), in 
upper right panel, we plot radial velocity Vr(r), in lower left 
panel, we plot A(r) and the in lower right panel, we plot 
Mach number Mx{r) [same as Mrir) in axial direction]. 
Matter with sound velocity a = 0.12272, injection veloc¬ 
ity Vin = 0.12319 and specific angular momentum A = 1.8 
is injected at Xin = 26.15. Each result has been obtained 
starting with an inviscid flow a = 0.0 and then increasing 
a up to 4.6 • 10“"^. In each case, standing shock solution 
was obtained. Note, that in this particular case, the shock 
disappears around a = 4.6 x 10~^ which is the critical a 
here. The shock location shifts outwards as described in COO 
and C96 with viscosity, and the angular momentum becomes 
more and more Keplerian. The results are indeed similar to 
Chakrabarti & Molteni (1995) solution of viscous isothermal 
flows. As viscosity is increased, the width of the shock also 
increases as it should be. For a different choice of injected 
parameters, the critical a could be very high as well (C96). 

In Fig. 2, we use the injected angular momentum to be 
same as the marginally stable value A = 1.8363 at the outer 
boundary located at Xin = 50. Injection velocity of matter 
and the sound speed are 0.157 and 0.0738 respectively. The 
stationary shock locations are seen in all the viscosity pa¬ 
rameter values. The critical viscosity in this particular case 
is much higher, and therefore the shock continues to persist 
for the range of a used. The angular momentum distribu¬ 
tion becomes closer to Keplerian much faster for higher val¬ 
ues of viscosity parameters. However, the flow remains sub- 
Keplerian throughout. By suitably adjusting the injection 
parameters at the outer boundary one obtains the solutions 
reaching a Keplerian disc. This requires very large number 
of smoothed particles. 

3. LOCATION OF STANDING SHOCKS IN THICK 
VISCOUS ACCRETION DISCS 

In this Section, we describe the formation of standing 
shocks in rotating, thick, axisymmetric, accretion discs onto 
a black hole. We introduce viscosity to the inviscid solutions 
described in detail by Molteni, Lanzafame & Chakrabarti, 
1994 (hereafter referred to as MLC94). In MLC94, the en¬ 


tire parameter space spanned by the specific energy e and 
specific angular momentum A was depicted which shows the 
boundary between shock and no-shock solutions in a flow 
in vertical equilibrium. It was observed that the theoreti¬ 
cal predictions based on flow models in vertical equilibrium 
regarding the shock locations are generally in good agree¬ 
ment with the simulation results when appropriate effects 
of turbulence is taken care of. 

As in the one-dimensional simulations presented in the 
previous section, we introduce viscosity gradually also in the 
two dimensional simulations, and in each case, we continue 
simulations till a stationary solution is reached. Secondly, 
unlike in MLC94, where the simulations were carried out 
with matter injected only in one quadrant, here we inject 
matter in two quadrants - both the upper and the lower side 
of the equatorial plane. As in MLC94, the accretion rate M 
being an eigenvalue of the problem, it is fixed by the choice 
of the input parameters (e. A). Therefore, matter density 
is chosen to be equal to unity. Since the accretion rate is 
an eigenvalue of the problem, it is automatically adjusted 
once the specific angular momentum and the specific energy 
are fixed at the outer boundary. This is true when a steady 
state is reached. When a steady state is not reached, such 
considerations do not apply. 

Figure 3 shows the XZ location of the SPH particles 
of the inviscid accretion disc model. Particles with specific 
angular momentum A = 1.6 and energy e = 0.001955219 
(which corresponds to injected radial velocity and the sound 
speed as Vin = 0.121122 and Uin = 0.058989 at Xi„ = 30 
respectively) are injected from the outer edge of the disc at 
A = 30. The results in the upper left panel is for inviscid 
flow. It is shown at t = 2900 when the flow has achieved 
a steady state. The total number of particles, whose size is 
h = 0.3, is 6146. A stationary shock is formed around X ~ 
5.2 and weaker oblique shocks are also seen (e.g., MLC94). 
At the two wings of the hot post-shock flow (basically inside 
the centrifugal barrier and the funnel wall, see, MRC96) a 
subsonic cooler outflow is present, symmetrically both up 
and down which carry away of the order of 5% of the total 
mass of the disc. No turbulent motions are evident in this 
model in the narrow post-shock subsonic flow. Subsequent 
to the shock, the hot flow becomes supersonic again before 
entering inside the black hole. 

In the other three panels we show the effect of the in¬ 
troduction of viscosity in the flow. Clearly, as in the one¬ 
dimensional case, higher viscosity causes higher differential 
angular momentum transport between the pre- and the post¬ 
shock solutions and as a result the shock is drifted away in 
the radial direction. The value of the viscosity parameter a 
is written on each panel. In the increasing order of viscosity, 
the shocks locations are Xs = 6.0, 7.5, and 10.5 respec¬ 
tively. Another point of interest: as viscosity is raised, the 
amount of outflowing matter in the wind is decreased. This 
is because of weakening of the centrifugal barrier of the in¬ 
going matter. Lower viscosity causes matter to bounce from 
the barrier and fly away as winds. In higher viscosity, higher 
turbulences are also seen to form in the post-shock flow. 
This is because more matter from higher elevation falls on 
the equatorial plane and convert their potential energy to 
turbulent energy. 

It is now widely believed (see, e.g., Chakrabarti, 1997 
and references therein) that the observed X-rays in black 
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hole candidates are produced and reprocessed in the cen¬ 
trifugal barrier supported regions very close to the black 
holes. The soft photons released from the pre-shock Keple- 
rian disc is reprocessed in the centrifugal pressure supported 
region (with or without shock) to become hard X-rays. Occa¬ 
sionally, radiations show quasi-periodic oscillations. Molteni, 
Sponholz & Chakrabarti (1996), and Ryu, Chakrabarti & 
Molteni (1997) suggest that these so-called QPOs observed 
in the compact objects could be due to the oscillations of the 
‘boundary layers’, i.e., the centrifugally supported denser re¬ 
gion. These oscillations form for those flow parameters where 
the theoretical analysis does not predict the formation of 
steady shocks, even though two sonic points exist. However, 
so far, no viscous flow was simulated. In Fig. 4a we show 
the first simulations of viscous flow shocks oscillate period¬ 
ically. In the upper-left, upper-right and lower-left panels 
roughly half of the cycle is shown where the shock location 
decreases monotonically. In the lower-right panel, the shock 
drifted again outward. The parameters are same as in Fig. 3, 
but the viscosity is higher: a = 2 x 10“®. The time period is 
around Tgpo = 4000 in units of 2GM/c^. Thus, in a stellar 
black hole of mass IOMq, the period would be around 0.4s, 
whereas in a supermassive black hole of mass lO’^M©, the 
period would be, 10®s or a few days. These time scales are 
comparable with the time scale of observed QPOs in black 
hole candidates. With a different disc input parameters (e.g., 
for different accretion rate or viscosity of the original disc) 
the periods will vary. Clearer evidence of the oscillation of 
the shock is shown in Fig. 4b, where the matter close to the 
equatorial place is collected and their Mach numbers are 
plotted. The behaviour of the shock on the equatorial plane 
is similar to what is seen in the one-dimensional simulations 
of Molteni, Sponholz & Chakrabarti, 1996. 

Fig. 4b: Mach number of the particles close to the equatorial 
plane is plotted against radial distance at different phases of os¬ 
cillation for the case shown in Fig. 4a 

Figure 5 shows the variation of the number of simula¬ 
tion particles as a function of time in the case shown in Fig. 
4a. Here the total number N and the number of sub-sonic 
particles Naub (presumably, participating in the Compton 
reprocessing of the soft-photons from the Keplerian disc) are 
shown. The amplitude modulation is signihcant: 50 percent 
variation in total particle number, and more than 700 per¬ 
cent variation in the sub-sonic particle number. We believe 
that this is important: the observed signihcant (10 — 100 per¬ 
cent) variations in QPO cannot be explained away by any 
means other than such a dynamical variation of the X-ray 
emitting region. 

4. CONCLUDING REMARKS 

In the present paper, we have numerically studied the 
behaviour of viscous, transonic flows. There are several im¬ 
portant conclusions: We show that the location of the cen¬ 
trifugally driven standing shock wave drifts away from the 
black hole and ultimately the shock disappears. This is valid 
both for one and two dimensional axisymmetric discs origi¬ 
nating from Keplerian flows. This result generalizes the ear¬ 
lier study of Chakrabarti & Molteni (1995) for isothermal 
hows done in one spatial dimensions. Second most impor¬ 
tant conclusion is that the viscous hows also show shock 
oscillations in regions of the parameter space where the 
steady shock is not possible even when two saddle type sonic 
points are present. This generalizes earlier work of Ryu, 


Chakrabarti & Molteni (1997) where non-viscous how was 
studied in two dimensions, and that of Molteni, Sponholz 
& Chakrabarti (1996) where non-viscous how was studied 
in presence of bremsstrahlung cooling. We believe that both 
these hndings are very important in explaining the observed 
X-ray radiations from the black hole candidates. The drift¬ 
ing shock solutions present the insight of how a Keplerian 
disc is actually formed out of an original sub-Keplerian how 
and how the centrifugal pressure supported region in the 
sub-Keplerian region around the black hole may act as the 
so-called Compton cloud. The dynamical oscillation of the 
centrifugal pressure supported shock wave (which, for the 
hrst time is shown here to be present even in a viscous tran¬ 
sonic disc) produces right frequency and the amplitude of 
modulation to be taken a serious explanation of the quasi- 
periodic oscillations. 
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Figure Captions 

Fig. 1 : Variation of density (upper-left), radial velocity (upper- 
right), angular momentum (lower-left) and Mach number (lower- 
right) with distance in a viscous transonic flow. Different curves 
are drawn for different viscosity parameters. As viscosity is in¬ 
creased, the shock is drifted outwards, until it disappears alto¬ 
gether beyond a critical viscosity parameter. The specific angular 
momentum with which matter is injected at the outer boundary 
is A = 1.8 

Fig. 2 : Same as in Fig. 1 except that the specific angular momen¬ 
tum of the injected flow is chosen to be the same as the marginally 
stable value of 1.8363. 

Fig. 3: Drifting of the steady shock location in a two dimensional, 
axisymmetric accretion flow as viscosity parameter is increased 
towards the critical value. See text for injected parameters. 

Fig. 4a: Oscillation of the shocks with time when a steady shock 
is not predicted, yet two sonic points exist in a viscous transonic 
flow. These solutions provide viable explaination to the quasi- 
periodic modulation of the X-rays in blak hole candidates. 

Fig. 5 : Oscillation of the total amount of matter in the accretion 
disc with a shock wave. Total particle number N and subsonic 
particle number Nguh plotted as functions of time. Parameters 
are same as those used to draw Fig. 4(a-b). 
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